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Abstract 

The problem of predicting the upper cut-off of the cosmic ray spec- 
trum is re-examined. Some aspects of extremely high energy interac- 
tions and their implications for the interpretation of giant air showers 
are also discussed. 



To my mother 



1 The Why of this Thesis 

The nature and origin of cosmic radiation has presented a challenge since 
its detection, some 80 years ago, in a series of pioneering ballon flights by 
Victor Hess [|l|. Already in 1938, it was clear from the observation of ex- 
tensive air showers (produced by energetic cosmic rays (CRs) interacting in 
the high atmosphere) @ that the CR-spectrum reached at least 10 16 eV, and 
continuously running monitoring indicated that the departure from isotropy 
was no greater than the statistical uncertainty of the measurements. Af- 
ter the discovery of the microwave background radiation (MBR) ||, it be- 
came obvious that ultra high energy particles undergo reactions with the 
relic photons (which they see as higly blue shifted), and so extremely high 
energy cosmic rays (EHECRs) cannot come from cosmologically large dis- 
tances |§].[] Assuming a cosmologically homogeneus population of sources, 
usually referred to as the universal hypothesis, this interaction produces the 
Greisen-Zatsepin-Kuz'min (GZK) cut-off around 5 x 10 19 eV. However, inge- 
nious installations with large effective areas and long exposure times needed 
to overcome the increasingly low flux, have raised the maximum observed 
primary particle energy to higher than 10 20 eV (see, e.g. , || for a recent 
review). In particular, the Fly's Eye experiment recorded the world's high- 
est energy cosmic ray event to date, with an energy proximately to 300 EeV 
M. The evidence for anysotropy of EHECR arrival directions is suggestive 
but statistically very weak. The analyses of the data sets of both Haverah 
Park J7| and AGASA |8j have suggested that the events with primary en- 
ergy > 40 EeV reached the Earth preferentially from the general direction of 
the Supergalactic Plane, a swath in the sky along which radio galaxies are 
clustered. However, the Fly's Eye group reported a small anisotropy to- 
wards the Galactic Plane at energies around 10 18 eV, confirmed later by the 

1 Usual CR- folklore makes the following distinction: ultra-high energy CRs are those 
above 10 18 eV while EHECRs refers to the ones with energies above 10 20 eV. History does 
not always choose the best variables, and I shall undistinguishably refer to them, specifing 
explicitly when speaking of events above 10 20 eV. 
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the AGASA experiment [|K]] (two analyses that did not reveal a significant 
correlation with the Supergalactic Plane). Recent observations dramatically 
confirms that the CR-spectrum does not end with the GZK cut-off Be- 
sides, the newest data world sample above 10 20 eV has no imprint of possible 
correlations with the Galactic Plane or Supergalactic Plane ||12||.p| 

Theoretical subtleties surrounding the production of particles with en- 
ergies below 10 20 eV have remained a puzzle. Acceleration in astrophysi- 
cal settings occurs when charged primaries, namely protons and heavy nu- 
clei, achieve high energies through repeated encounters with moving magne- 
tized plasmas ("bottom up" mechanisms) [Q. Preferred sites are supernova 
shocks, galactic wind termination shocks, or relativistic shocks associated 
with active galactic nuclei (AGNs) and radio galaxies (known to be powerful 
sources of radio, gamma radiation, and super-relativistic jets). Although in 
this last case, acceleration up to around 10 21 eV seems to be possible by 
stretching the reasonable values for the shock size and the magnetic field 
strength at the shock somewhat, for the highest events, there seem to be no 
suitable extragalactic objects such us AGNs or active galaxy clusters near 
the observed arrival direction and within a maximum distance of about 50 
Mpc (the estimated angular deflection for a proton of 10 20 eV is 2.1°). The 
gamma ray burst (GRB) phenomenon, which itself is also an oustanding 
mistery in astrophysics, has been considered as one of the candidates for 
the origin of EHECRs |15[ . GRBs are generally related with catastrophic 



(explosive) events during which CR particles may be shocked accelerated to 
extreme energies, besides, the total energy and ocurrence rate of GRBs are in 
rough agreement with the production rate of EHECRs. However, the model 



is not free of problems, as can be see, for instance, in [|16|j . 

The difficulties encountered by conventional acceleration mechanisms in 
accelerating particles to the highest observed energies have motivated sug- 
gestions in the sense that the underlying production mechanism could be, 
perhaps, of non-aceleration nature. In the so-called "top down" models, 
charged and neutral primaries are produced at extremely high energies, typ- 
ically by quantum mechanical decay of supermassive elementary X particles 
related to grand unify theories |L7]]. These exotic particles have also been 
suggested as primaries 



2 A possible correlation between compact radio quasars and the five most energetic CRs 
has been already reported [Of . 
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An alternative (perhaps mischievously) explanation for the EHECRs re- 
quires the breakdown of local Lorentz invariance. Namely, tiny departures 
from Lorentz invariance, too small to have been detected otherwise, might 
affect elementary particle kinematics in such a way that some hadronic reso- 
nances which are inestable at low energies would become stable at very high 
energies. Therefore the GZK cut-off can be relaxed or removed |19[ . It could 
be found some even more "exotic" : The recently detected energetic particles 
may have shortcut their journey traversing a wormholef] The question then 
could be not from where do the rays come, but from when!!! 

Deepening the mistery, the particle identity is still unknown. The Fly's 
Eye analysis suggested a transition from a spectrum dominated by heavy 
nuclei to a predominantly light composition above a few times 10 19 eV. How- 
ever, it may be worth mentioning that proton induced air shower Monte 
Carlo does not fit the shower development of the highest event with high 
precision. A primary heavy nucleus fits more closely the shower development 
and neutrinos cannot be excluded as atypical primaries neither The 
situation is not settled down because fluctuations in the shower development 
are known to be large. On the one hand, heavy nuclei have their own mer- 
its because they can be deflected considerably by the galactic magnetic field 
which relaxes the source direction requirements, and they can be accelerated 
to higher terminal energies because of their higher charge. On the other 
hand cosmic ray nuclei with energies in excess of 2 x 10 20 eV cannot origi- 
nate from sources beyond 10 Mpc |22|, |23| . Weakly interacting particles such 
as neutrinos will have no difficulty in propagating through the intergalactic 
medium, namely the corresponding mean free path I„si 4 x 10 28 cm is just 
above the present size of the horizon, Hq 1 m 10 28 cm pl| . Many sources of 
extremely high energy neutrinos are known, they can be produced in astro- 
physical objects by the decay of pions or kaons generated as subproducts of 
proton-photon interactions during the acceleration processes |25|] or else by 
"top down" mechanisms [I7|. In the latter, the neutrino flux might extend 



even up to 10 25 - 10 28 eV. However, at these energies the atmosphere is still 
transparent to neutrinos, and most of them will impact on Earth Cos- 



mic ray neutrino annihilation onto relic neutrinos in the galactic halo was 
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Naively, a wormhole is a tunnel in the topology of the space-time from where in-coming 



causal curves can pass through and become out-going on the other side 1 20 . Whether such 



wormholes arc actually allowed by the laws of physics is currently unknown. 
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suggested as an alternative source of EHECRs P7| , but as the previous ones, 
this is not a truly satisfactory solution . 

Summing up, the GZK cut-off provides an important constraint on the 
proximity of EHECR (nucleons and nuclei) sources, thus, the origin of the 
events with energies above 100 EeV has became one of the most pressing 
questions of cosmic ray astrophysics. 

This summary is ordered as follows. The next section is devoted to revisit 
the effect of MBR on EHECRs using the continuous energy-loss (CEL). It is 
a matter of fact, that even when transport equations may be as complicated 
as one can afford, some simple assumptions makes results agree, within a 
typical few percent, with those genereted by numerical code, thus allowing 
for the isolation of essential physics. This is exactly what is needed here 
in order to study the the footprint left by the relic photons in the energy 
spectrum and how can this be used to extract testeable predictions. In Sec. 
3, using polarization measurements and synchrotron emission as a trace, we 
analyze concrete candidates for testing the acceleration models that involve 
hot spots and strong shock waves as the source of the ultra-high energy com- 
ponent of the CR-spectrum. Sec. 4 will deal with the influence of different 
hadronic models on extensive air showers. The hadronic models considered 
are those implemented in the well-known QGSJET and SIBYLL pO 



event generators. The different approaches used in both codes to model the 
underlying physics are analyzed using computer simulations performed with 
the program AIRES . 



2 En route to us, here on Earth 

The aim of this section is to present a simplify model of the propagation of 
EHECRs which reproduces the essential elements of the highly sophisticated 
computer calculations to within a few percent. 

2.1 Energy attenuation length 
2.1.1 Nucleons 

If one assumes that the highest energy cosmic rays are indeed nucleons, 
the fractional energy loss due to interactions with the cosmic background 
radiation at temperature T and redshift z = 0, is determined by the integral 
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of the nucleoli loss energy per collision times the probability per unit time 
for a nucleon collision moving through an isotropic gas of photons . This 
integral can be explicitly written as follows, 

IdE c ^ f» . v T m , "H m 

where iy r is the photon energy in the rest frame of the nucleon, and Kj is 
the average fraction of the energy lost by the photon to the nucleon (in the 
jth reaction channel) in the laboratory frame (i.e. the frame in which the 
microwave background radiation is at « 3K). The sum is carried out over 
all channels, n(w)dw stands for the number density of photons with energy 
between w and dw, following a Planckian distribution J33| at temperature 
T, <jj(w r ) is the total cross section of the (jth) interaction channel, T is the 
usual Lorentz factor of the nucleon, and w m is the maximum energy of the 
photon in the photon gas. 

Thus, the fractional energy loss is given by 



1 dE _ ckT 
E~dT " ~27r 2 r 2 (cft) 3 



dw r aj(w r ) Kj w r ln(l - e -^/2VhT^ ^ 



where k and h are Boltzmann's and Planck's constants respectively, and Wq. 
is the threshold energy for the jth reaction in the rest frame of the nucleon. 

The characteristic time for the energy loss due to pair production at 
E > 10 19 eV is t ~ 5 x 10 9 yr p4| and therefore it does not affect the spectrum 
of nucleons arriving from nearby sources. Consequently, for nucleons with 
E > 3 x 10 19 eV (taking into account the interaction with the tail of the 
Planck distribution), meson photoproduction is the dominant mechanism for 
energy loss. Notice that we do not distinguish between neutrons and protons; 
in addition, the inelasticity due to the neutron (3 decay is negligible. 

In order to determine the effect of meson photoproduction on the spec- 
trum of cosmic rays, we first examine the kinematics of photon-nucleon in- 
teractions. Assuming that reactions mediated by baryon resonances have 
spherically symmetric decay angular distributions []3~5[| , the average energy 



loss of the nucleon after n resonant collisions is given by 

K(m Ro ) = l--Y[(l+ \ 2 M \ (3) 



ni 



R. 



5 



where ttlr denotes the mass of the resonant system of the chain, m,M the 
mass of the associated meson, = yfs is the total energy of the reaction 
in the centre of mass, and m.p the mass of the nucleon. It is well established 
from experiments |36| that, at very high energies (yfs above ~ 3 GeV), 
the incident nucleons lose one-half their energy via pion photoproduction 
independently of the number of pions produced (leading particle effect). 

In the region dominated by baryon resonances, the cross section is de- 
scribed by a sum of Breit-Wigner distributions (constructed from the exper- 
imental data in the Table of Particle Properties |37|]) considering the main 
resonances produced in N^y collisions with ttN, ttttN, and ttKN final states. 
For the cross section at high energies we used the fits to the high-energy cross 
section OtotaiM made by the CERN, DESY HERA, and COMPAS Groups 



3jJ. In this energy range, the cr to tai(^7) is to a good approximation identical 



to <7total(P7)- 

The numerical integration of Eq. ([|) is performed taking into account 
the aforementioned resonance decays and the production of multipion final 
states at higher centre of mass energies (y/s ~ 3 GeV). A \ 2 fit of the numer- 
ical results of equation (0) with an exponential behaviour, Aexp{— B/E}, 
proposed by Berezinsky and Grigor'eva |38| for the region of resonances, 
gives: A = (3.66 ± 0.08) x 10- 8 yr _ \ B = (2.87 ± 0.03) x 10 20 eV with a 
X 2 /dof = 3.9/10. At high energies the fractional energy loss was fitted with 
a constant, C = (2.42 ±0.03) x 10~ 8 yr _1 |39 |. These results differ from those 



obtained in , due to a refined expression for the total cross section. From 
the values determined for the fractional energy loss, it is straightforward to 
compute the energy degradation of EHECRs in terms of their flight distance. 

In Fig. 1 we compare the energy attenuation length of nucleons obtained 
with the CEL approximation and the Monte Carlo simulation taken from Ref. 
pO|j . Although the energy loss is dominated by single pion photoproduction 



hj. |4T|1 , one can see that at higher energies the multiple production of pions 
becomes quite important. Notice that, independently of the initial energy of 
the nucleon, the mean energy values approach to 100 EeV after a distance 
of « 100 Mpc. 



2.1.2 Nuclei 

The relevant mechanisms for energy losses that extremely high energy nu- 
clei suffer during their trip to the Earth are photodisintegration and hadron 
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Figure 1: Energy attenuation length of nucleons in the intergalactic medium. 
Solid line stands for the CEL approximation, while dots for Monte Carlo 
simulation which considers just single pion production ||40|| . 



photoproduction (which has a threshold energy of ~ 145 MeV, equivalent to 
a Lorentz factor of 10 11 , above the range treated in this Thesis). 

Following the conventions of Eq. (H), the disintegration rate of a nucleus 
of mass A with the subsequent production of i nucleons is given by the 
expression 



42 



RAi 



2T 2 



00 , n[w) 
dw — — 



2Yw 



dw r W r <JAi(w r 



(4) 



with OAi the cross section for the interaction. Using the expressions for the 



cross section fitted by Puget et al. |E| , it is possible to work out an analytical 
solution for the nuclear disintegration rates ||44|| . After adding over all the 
possible channels for a given number of nucleons, one obtains the effective 
nucleon loss rate. The effective 56 Fe nucleon loss rate obtained after carrying 
out these straightforward but rather lengthy steps can be parametrized by, 



R(T) = 3.25 x l(r 6 r- a643 exp(-2.15 x 10 10 /r) s' 1 



(5) 
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if T G [1. x 10 9 ,3.68 x 10 10 ], and 

i?(r) = i.59xio- 12 r-°- 0698 s- 1 (6) 

if T G [3.68 x 10 10 , 1. x 10 11 ]. It is noteworthy that knowledge of the iron 
effective nucleons loss rate alone is enough to obtain the corresponding value 



of R for any other nuclei [43 



The emission of nucleons is isotropic in the rest frame of the nucleus, 
and so the averaged fractional energy loss results equal the fractional loss in 
mass number of the nucleus, viz., the Lorentz factor is conserved. Because 
of the position of 56 Fe on the binding energy curve, it is considered to be a 
significant end product of stellar evolution, and higher mass nuclei are found 
to be much rarer in the cosmic radiation. Thus, hereafter, when speaking 
of sources of heavy nuclei we shall be thinking on iron nuclei. The relation 
which determines the attenuation length for its energy is then, 

E = E g e- R ^ t / 56 (7) 

where E g denotes the energy with which the nuclei were emitted from the 
source, and Y = E g /56. 

In Figs. 2 and 3 we plot the final mass (i.e. A) and energy E of the 
heaviest surviving fragment as a function of the distance for initial iron nuclei 
(a). For comparison, it is also super- imposed a Monte Carlo simulation which 
include the rates just discussed and also the pair creation processes (b) [pfl . 
One can see that nuclei with Lorentz factors above 10 10 cannot survive for 
more than 10 Mpc (recall thet 10 Mpc correspond to ~ 10 15 s), and for these 
distances, the CEL approximation always lies in the region which includes 
95% of the simulations. 

At To < 5 x 10 9 the pair creation losses start to be relevant, signifi- 
cantly reducing the value of the Lorentz factor as the nucleus propagates 
distances of (9(100 Mpc). The effect has a maximum for T ~ 4 x 10 9 but 
becomes small again for T < 10 9 , for which appreciable effects only appear 
for cosmological distances (> 1000 Mpc), see for instance PB|. The efect of 



neglecting pair creation losses translates into keeping T = r constant during 
the propagation, and this enhances the photodisintegration rates and then 
reduces < A > more rapidly. The divergences between (b) and (c) are less 
pronounced in Fig. 3. Namely, for T = 4 x 10 9 the average energies with and 
without pair creation proccess (obtained with the Monte Carlo simulations) 
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Figure 2: Evolution of the mass number A of the heaviest fragment surviving 
photodisintegration vs. travel time for different injection energies (Tq = 
4 x 10 9 , To = 2 x 10 10 ). (a) stands for CEL approximation. It is also included 
< A > obtained from Monte Carlo simulations (with (b) and without (c) pair 
creation processes) for comparison. The region between the two dashed lines 
includes 95% of the simulations. This gives a clear idea of the range of values 
which can result from fluctuations from the average behaviour. 



are similar up to ~ t = 10 16 s while the < A > values sizeably differ for 
t = 3 x 10 15 s onwards. When comparing with the CEL the differences arise 
even earlier (« 70 Mpc and 30 Mpc respectively). This compensation is due 
to a partial cancellation between the effects of the evolution of F and of A 
in the values of the final energy (E = m p T A) , since neglecting pair creation 
losses does not allow F to decrease but make instead A to drop faster. 

In Fig. 4 we have shown the relation between the injection energy and the 
energy at a time t for different propagation distances. The graph indicates 
that the final energy of the nucleus is not a monotonic function. It has a 
maximum at a critical energy and then decreases to a minimum before rising 
again as T rises, as was first pointed out by Puget et al. [43|. The fact that 
the energy E is a multivaluated function of E g leads to a pile-up in the energy 
spectrum. Moreover, this behaviour enhances a hidden feature of the energy 
spectrum for sources located beyond 2.6 Mpc: A depression that preceeds a 
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Figure 3: Energy of the surviving nuclei (T = 4 x 10 9 , T = 2 x 10 10 ) 
vs. propagation distance obtained with the CEL (a). It is also included the 
energy attenuation length computed with Monte Carlo simulations, while 
(b) include both photodisintegration and pair creation processes, (c) just the 
losses due to photodisintegration. Again, the region between dashed lines 
includes the 95 % of the simulations. 

bump that would make the events at the end of the spectrum (just before 
the cut-off) around 50% more probable than those in the depressed region. 

2.2 Modification of the cosmic ray spectrum 
2.2.1 A hump follow by a cut-off 

Let us begin with the modification that the MBR produces in the EHECR 
(nucleon) spectrum. The evolution of the spectrum is governed by the bal- 
ance equation 

ON d(b(E)N) „ _ 2Ar ^ 

that takes into account the conservation of the total number of nucleons 
in the spectrum. In the first term on the right, b(E) is the mean rate at 
which particles lose energy. The second term, the diffusion in the MBR, 
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Figure 4: Relation between the injection energy of an iron nucleus and the 
final energy of the photodisintegrated nucleus for different values of the prop- 
agation distance (from grey to black, 20 Mpc, 10 Mpc, 3.5 Mpc, and 3 Mpc.) 



is found to be extremely small due to the low density of relic photons and 
the fact that the average cosmic magnetic field is less than 10~ 9 G EB| and 
is neglected in the following. The third term corresponds to the particle 
injection rate into the intergalactic medium from some hypothetical source. 
Since the origin of cosmic rays is still unknown, we consider three possible 
models: 1) the universal hypothesis, which assumes that cosmic rays come 
from no well-defined source, but rather are produced uniformly throughout 
space, 2) single point sources of cosmic rays, and 3) sources of finite size 
approximating clusters of galaxies. In all the cases it has been assumed that 
nucleons propagate in a straight line through the intergalactic medium due 
to the reasons mentioned above that allow us to neglect the diffusion term. 

There exists evidence that the source spectrum of cosmic rays has a power- 
law dependence Q(E) = kE' 1 (see for instance ||7|])- With the hypothesis 
that cosmic rays are produced from sources located uniformly in space, with 
this power law energy dependence (which implies a steady state process), the 
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solution of Eq. (||) is found to be 

, s Q(E) E 

For the case of a single point source, the solution of equation (§) reads 



N{E ' t) = biE) Je Q&vVdEv (10) 



b(E) 

with 

rE a dE 



t' = t- / (11) 
Je b(E) V ; 

and Eg the energy of the nucleon when emitted by the source. The injection 
spectrum of a single source located at to from the observer can be written 
as Q(E,t) = kE~^ 8(t — to) and for simplicity, we consider the distance as 
measured from the source, that means to — 0. At very high energies, i.e. 
where b(E) = C E (with C the constant defined above in the discussion of 
fractional energy loss) the total number of particles at a given distance from 
the source is given by 

w <*«>~iCgjjCV'(«4 h §H (12) 

or equivalently, 

N(E,t) ^ KE~^e-^- 1)ct . (13) 

This means that the spectrum is uniformly damped by a factor depending 
on the proximity of the source. 

At low energies, in the region dominated by baryon resonances, the 
parametrization of b(E) does not allow a complete analytical solution. How- 
ever using the change of variables, 

~ f E <> dE 

t= / 14 
Je b(E) 

with Eg = £(E,t) and dt = dE g /b(E g ), we easily obtain, 

N(E, t) = —j o £(E, t)-' 5(t - t) b[£(E, t)} dt, (15) 
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and then, the compact form, 



N(E,t) = -^E^b{E g ). (16) 

In our case, where we have assumed an exponential behaviour of the fractional 
loss energy, E g is fixed by the constraint: At — Ei (B/E) + Ei (B/E g ) = 0, 
Ei being the exponential integral ||49|| , and B the constant defined above in 



the parametrization of Berezinsky and Grigor'eva. 

Studies are underway of the case of a nearby source (i.e. within about 
100 Mpc) for which the fractional energy loss AE/E is small. Writing 
E g = E + AE, and neglecting higher orders in AE, it is possible to obtain 
an analytical solution and the most relevant characteristics of the modified 
spectrum. Therefore, an expansion of b(E + AE) in powers of AE in Eq. 
( |TTD allows one to obtain, an expression for AE, 



AE w < exp 



tb(E) (E + B) 
E~ 2 



E2 d7) 



{E + B)' 



To describe the modification of the spectrum, it is convenient to introduce 
the factor rj, the ratio between the modified spectrum and the unmodified 
one, that results in 

f E + AE\^ b(E + AE) , . 

'=(-6-) V' (18) 

Equations flT7|) and ( |T8|) describe the spectral modification factor up to en- 
ergies of ~ 95 EeV (~ 85 EeV) with a precision of 4% (9%) for a source 
situated at 50 Mpc (100 Mpc). 

In Fig. 5 we plot the modification factors for the case of near sources with 
power law injection (7 = 2.5) to compare with the corresponding results of 



3q| . The most significant features of rj are the bump and the cut-off. It has 
been noted in the literature that the continuous energy loss approximation 
tends to overestimate the size of the bump |50], [51]]. Fig. 5 shows how the 
bump is less pronounced in our treatment. As the bump is a consequence 
of the sharp (exponential) dependence of the free path of the nucleon on 
energy, we attribute these differences to the replacement of a cross section 
approximated by the values at threshold energy PS] , by a more detailed 



13 



5 3.5 



O 



c 
o 



3 - 



£ 2.5 



O 



1 .5 



0.5 - 





this work 




Berezinsky et ol 


50 Mpc 




this work 




- ,* E 


erezinsky et ol 


100 Mpc 


i iii 


i i \ 





10 



EteV] 



Figure 5: Modification factor of single-source energy spectra for different 
values of propagation distance and power law index 7 = 2.5. 



expression taking into account the main baryon resonances that turn out to 
be important. 

Another alternative for the particle injection comes from clusters of galax- 
ies. These are usually modelled by a set of point sources with spatially uni- 
form distribution, although in reality the distribution of galaxies inside the 
clusters is somewhat non-uniform. In our treatment we shall assume that the 
concentration of potential sources at the center of the cluster is higher than 
that in the periphery, and we adopt a spatial gaussian distribution. With 
this hypothesis, the particle injection rate into the intergalactic medium is 
given by 

Q(E, t) = * f° JjZ- 8(t - T) exp { ~ (T ~ to)2 | dT (19) 
J -00 V2 7T a [ 2a z J 

A delta function expansion around to, with derivatives denoted by lower case 
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Roman superscripts, 

8{t — T) = 5(t - t ) + (t - t ) (T - t ) + ^ (t — t )(T — t ) 2 + ... (20) 
leads to a convenient form for the injection spectrum, which is given by, 

Q(E, t') = kET* [S(f - t ) + ^ 8^(t' - to) + ^ <5^(t' - to) + • • •] (21) 

From Eqs. fliCf ) and (|TT]), it is straightforward to compute an expression for 
the modification factor, 



E-^HEa) { a 2 A 2 e- 2B l E <> a 4 A 4 e" 4B/ ^ i 



where 
and 



(22) 

F x {E g ) = 2B 2 E; 2 + (2 - ^BE' 1 + (1 - 7 ) 2 , 



F 2 (£ 5 ) = 245 4 E- 4 + (4- 50 7 )5 3 E- 3 + (35 7 2 - 257 + 8)5^; 



3 i on, ,2 ir„ , ^ \ n zt 1 — 1 , /i „,\4 



+ (-10 7 d + 20y - 157 + 4)5£g + (1 - 7 ) 4 . 

The modification factor for the case of extended sources described by 
Gaussian distributions of widths 2 and 6 Mpc at a distance of 18.3 Mpc 
is shown in Fig. 6. These may be taken as very crude models of the Virgo 



cluster |52j , assuming that there is no other significant energy loss mechanism 
for cosmic rays traversing parts of the cluster beyond those due to interactions 
with the cosmic background radiation. 

A power law injection ( 7 = 2.5) was used again, as in the pointlike case. 
The curves can be understood qualitatively. Both peaks are in about the 
same place, occurring around the onset of pion photoproduction. The broader 
source distribution reflects the losses suffered by the more remote part of the 
distribution in traversing a greater distance to us. 

The high energy cosmic ray spectrum as measured by Fly's Eye detector, 
both in stereo and monocular modes, has been fitted to Eq. (|I1J) and k E~ j x 
Eq. (|22|). The results are shown in Table I. The first two rows do not include 



the experimental energy resolution [E31 while the third one does 



4 It is important to stress that the if the most energetic event is included, the single 
source as well as the extended source model, cannot explain the data. 
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Figure 6: Modification factor of extended-source energy spectrum for a prop- 
agation distance ~ 18.3 Mpc and power law index 7 = 2.5. 

2.2.2 A depression before a bump on the EHECR nuclei spectrum 

The photodisintegration process results in the production of nucleons of ul- 
trahigh energies with the same Lorentz factor of the parent nucleus. As a 
consequence, the total number of particles is not conserved during propa- 
gation. However, the solution of the problem becomes quite simple if we 
separately treat both the evolution of the heaviest fragment and those frag- 
ments corresponding to nucleons emitted from the travelling nuclei. The 
evolution of the differential spectrum of the surviving fragments is again 
governed by a balance equation that takes into account the conservation of 
the total number of particles in the spectrum. Using the formalism presented 
in the previous section, and considering the case of a single source located 
at to from the observer, with injection spectrum Q(E g ,t) = nEg 1 5{t — to), 
the number of particles with energy E at time t is given by, 

kE~ 1+1 

N(E, t)dE = — 9 - — dE, (23) 
E 

with Eg fixed by the constraint ([?]). 
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Table 1: Results of the fits performed by Dova et al. [53,54]. 



source 


7 


log K 


distance [Mpc] 


a [Mpc] 


single 


3.27 ± 0.03 


29.64 ± 0.02 


11312 




extended 


3.27 ± 0.03 


29.64 ± 0.02 


93 ± 15 


30 ± 10 


extended 


3.27 ± 0.02 


29.30 ± 0.45 


119 ± 33 


35 



Let us now consider the evolution of nucleons generated by decays of 
nuclei during their propagation. For Lorentz factors less than 10 11 and dis- 
tances less than 100 Mpc the energy with which the secondary nucleons are 
produced is approximately equal to the energy with which they are detected 
here on Earth. The number of nucleons with energy E at time t can be ap- 
proximated by the product of the number of nucleons generated per nucleus 
and the number of nuclei emitted. When the nucleons are emitted with en- 
ergies above 100 EeV the losses by meson photoproduction start to become 
significant. However, these nucleons come from heavy nuclei with Lorentz 
factors T > 10 11 which are completely disintegrated in distances less than 
10 Mpc. Given that the mean free path of the nucleons is about A„ ~ 10 
Mpc, it is reasonable to define a characteristic time r r given by the moment 
in which the number of nucleons is reduced to 1/e of its initial value A . In 
order to determine the modifications of the spectrum due to the losses which 
the nucleons suffer due to interactions with the relic photons, we assume that 
the iron nucleus emitted at t = t is a travelling source which at the end of 
a time r r has emitted the 56 nucleons altogether. In this way the injection 
spectrum of nucleons (T « 10 11 ) can be approximated by, 

q{E G ,t) = KA^ +l E G H(t-T T ), (24) 

where A Q is the mass of the initial nucleus and the energy with which the 
nucleons is generated is given by Eq = E g /A . 

The number of nucleons with energy E at time t is given by, 

n (E, t)dE = dE (25) 

E 

and the relation between injection energy and the energy at time t remains 
fixed by the relation, A(t — r r ) — Ei (B/E) + Ei (B/E G ) = 0, Ei being the 
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Figure 7: Modification factors of heavy nuclei sources located at 20 Mpc. 

exponential integral, and A, B the parameters of the fractional energy loss 
of nucleons fitted in the previous section. 

In Fig. 7 it is shown the modification factors for the case of sources of iron 
nuclei (propagation distance 20 Mpc) together with the spectra of secondary 
nucleons. It is clear that the spectrum of secondary nucleons around the 
pile-up is at least one order of magnitude less than the one of the surviving 
fragments. In Figure 8 we have plotted the modification factor for a source of 
nuclei located at 2.8 Mpc. It displays a bump and a cut-off and, in addition, 
a depression before the bump. A detailed analysis of the variation of these 
features with the propagation distance and the spectral index was presented 
It is important to stress that the mechanism that produces the pile- 



in 



up, which can be seen in Figs. 7 and 8, is completely different to the one 
that produces the bump in the case of nucleons (Figs. 5 and 6). 

In this last case, the photomeson production involves the creation of new 
particles that carry off energy, yielding nucleons with energies ever closer to 
the photomeson production threshold. This mechanism, modulated by the 
fractional energy loss, is responsible for the bump in the spectrum. The cut- 
off is a consequence of the conservation of the number of particles together 
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Figure 8: Modification factor of single source energy spectrum for a propaga- 
tion distance of 2.8 Mpc, assuming a differential power law injection spectrum 
with spectral index 7 = 2.5. 

with the properties of the injection spectrum E~^dE g < 00). 

th " 

In the case of nuclei, since the Lorentz factor is conserved, the surviving 
fragments see the photons of the thermal background always at the same 
energy. Then, despite the fact that nuclei injected with energies over the 
photodisintegration threshold lose energy by losing mass, they never reach 
the threshold. The observed pile-up in the modification factors is due solely 
to the multivalued nature of the energy at time t as a function of the injec- 
tion energy: Nuclei injected with different energies can arrive with the same 
energy but with different masses. 

It is clear that, except in the region of the pile-up, the modification fac- 
tor i] is less than unity, since rj = (E/E g )' y ~ l . This assertion seems to be in 
contradiction with the conservation of particle number. Actually, the con- 
servation of the Lorentz factor implies, 

KEpdE g \ T = N(E,t)dE\ r (26) 

in accord with the conservation of the number of particles in the spectrum. 
Moreover, the condition (|2"6"D completely determines the evolution of the en- 
ergy spectrum of the surving fragments (p3|). Note that in order to compare 
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the modified and unmodified spectra, with regard to conservation of particle 
number, one has to take into account that the corresponding energies are 
shifted. As follows from fl2"S|), the conservation of the number of particles in 
the spectrum is given by, 

/ N{E,t)dE= kET 1 dE g (27) 

with E th and E nth the threshold energies for photodisintegration, and photo- 
pion production processes respectively. 

Let us now return to the analysis of Fig. 4 in relation to the depres- 
sion in the spectrum. In the case of a nearby iron source, located around 3 
Mpc, and for injection energies below the multivalued region of the function 
E(E g ), E is clearly less than E g and, as a consequence the depression in the 
modification factor is apparent. Then, despite the violence of the photodis- 
integration process via the giant dipole resonance, for nearby sources none of 
the injected nuclei are completely disintegrated yielding this unusual depres- 
sion before the bump. For a flight distance of 3 Mpc, the composition of the 
arrival nuclei changes from A = 50 (for T « 10 9 ) to A = 13 (for F « 10 11 ). 
However, the most important variation takes place in the region of the bump, 
where A runs from 48 to 13, being heavy nuclei of A = 33 the most abun- 
dant. For propagation distances greater than 10 Mpc one would expect just 
nucleons to arrive for injection energies above 9 x 10 20 eV. In this case the 
function becomes multivalued below the photodisintegration threshold and 
then there is no depression at all. For an iron source located at 3.5 Mpc, the 
depression in the spectrum is almost invisible ((9(1%)), in good agreement 
with the results previously obtained by Elbert and Sommers using Monte 
Carlo simulation f26fl . 

The structure of the spectrum, recently published by AGASA experiment 
(Fig. 9) fil l, seems to suggest the presence of a bump around 4 x 10 19 eV; the 



central points may be unvailing a second bump around 2 x 10 20 eV. Whatever 
the source of the highest energy CRs, because of their interactions with the 
2.7K relic radiation and magnetic fields in the universe, those reaching the 
Earth will have spectra, composition and arrival directions affected by prop- 
agation. Within the universal hypothesis, nucleons with energies above 100 
EeV are pushed below this limit, and as a consequence, particles which origi- 
nally had a higher energy pile-up forming a bump around 5 x 10 19 eV followed 
by a sharped cut-off. While the origin of the highest energy cosmic rays is 
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Figure 9: Energy spectrum of EHECRs measured by AGASA experiment. 
The spectrum is normalized to the first event (E = 10 18 - 55 eV). 

still uncertain, it is not necessary to invoke top down mechanisms to explain 
the existing data. It seems likely that shock acceleration at Fanaroff-Riley 
radio galaxies, and Starburst galaxies located around 3 Mpc away (recall that 
nucleus energy, for T > 10 11 is degraded even faster than that of protons) 
can account for the existing data at the end of the CR-spectrum (energies 
above 100 EeV). The existing experimental data are insufficient to draw any 
definitive conclusions, and we are obliged to present this idea as a hypothesis 
to be tested by experiment. 

3 The southern CR-sky 

Radio galaxies considered as sources of CR beyond the GZK cut-off must 
be quite close to our own galaxy (z < 0.03). There are just a few objects 
within this range and, consequently (if these are responsable of the end of 
the CR-spectrum), the expected EHECR distribution on the sky should be 
highly anisotropic. Assuming that the intergalactic magnetic field is near 
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1 nG, it should be expected an excess of CR detections at energies larger 
than 5 x 10 19 eV within regions of angular radius 9 < 15° and centered 
at the positions of the nearest active galaxies. According to this picture, 
the southern CR-sky should be dominated by three outstanding sources: 
Centaurus A (the nearest radio galaxy) which may provide the most energetic 
particles detectable on Earth, Pictor A (a strong source with a flat radio 
spectrum) which would contribute with the larger CR flux p5| , and PKS 
1333-33 which might be a source of events similar to those recently detected 
in the Northern Hemisphere. In addition to these sources, there are other 
two southern candidates, Fornax A (z = 0.057) and PKS 2152-69 (z = 
0.027), which could provide contributions to the CR flux above the cut-off.0 
An approximate theoretical picture of the up to now unexplored ultra high 
energy CR southern sky is consequently at our disposal. In this section 
we shall discuss the main features of the radio-galaxies Cen A and PKS 
1333-33. Before doing so, let us briefly summarize diffuse shock acceleration 
mechanisms. 



3.1 Bottom up models 

The hot spots of extended radio sources are regions of strong synchrotron 



emission 



56| . These regions are produced when the bulk kinetic energy of 
the jets ejected by a central active source (supermassive black hole + ac- 
cretion disk) is reconverted into relativistic particles and turbulent fields at 
a "working surface" in the head of the jets |57|]. The speed with which 
the head of a jet advances into the intergalactic medium of particle density 
n e can be obtained by balancing the momentum flux in the jet against the 
momentum flux of the surrounding medium. Measured in the frame comov- 
ing with the advancing head, t>h ~ Vj [1 + (ne/rij) 1 / 2 ] -1 , where rij and Vj are 
the particle density and the velocity of the jet flow, respectively. Clearly, 
fj > t>h for n e > rij, in such a way that the jet will decelerate. The result is 
the formation of a strong collisionless shock, which is responsible for particle 
reacceleration and magnetic field amplification. The acceleration of parti- 
cles up to ultrarelativistic energies in the hot spots is the result of repeated 
scattering back and forth across the shock front. The particle deflection in 
this mechanism is produced by Alfven waves in the turbulent magnetic field. 

5 Concerning nuclei sources, it should be mentioned NGC 253. 
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This process has been studied in detail by Biermann and Strittmatter . 
Assuming that the energy density per unit of wave number of MHD turbu- 
lence is of Kolmogorov type, i.e. I(k) a k~ s with s = 5/3, the acceleration 
time scale for protons given by: 

\~l u (A) (28) 

where /3 JET is the jet velocity in units of c, u is the ratio of turbulent to 
ambient magnetic energy density in the hot spot (of radius R HS ), and B is 
the total magnetic field strength. The acceleration process will be efficient 
as long as the energy losses by synchrotron radiation and photon-proton 
interactions do not become dominant. Considering an average cross section 
<7 7P for the three dominant pion-producing interactions |59[| , 7 + p — > p + 7r°, 
7 + p — > n + iT + , 7 + p— > p + tt + + if- , we get the time scale of the energy 
losses within a certainty of 80%: 




6 7r mi c 



a T mlB 2 (l + a A 



3 

E- 1 (29) 



where a stands for the ratio of photon to magnetic energy densities, a T is 
the classical Thomson cross section, and A gives a measure of the relative 
strength of 7p interactions against the synchrotron emission. Biermann and 



Strittmatter |58| have estimated A f» 200, almost independently of the source 
parameters. The most energetic protons injected in the intergalactic medium 
will have an energy that can be obtained by balancing the energy gains and 
losses: 

£ max = 7.8 x 10 5 R^' 2 BZT (1 + A a)" 3 / 4 EeV, (30) 

here -B_5 is the magnetic field in units of 10~ 5 G. 

The diffusive shock acceleration process in the working surface leads to a 
power law particle spectrum [BD|, 



Q(E) = k ETi (E <E< £ max ) (31) 

where k = (no/7 ~ 1) ^0 7 f° r -^max ^> E , and n is the particle density in 
the source. 
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3.2 Centaurus A 



The most attractive features of Centaurus A (a complex and extremely pow- 
erful radio source identified at optical frequencies with the galaxy NGC 5128) 
as a possible source of CRs are its active nature, its large, non-thermal ra- 
dio lobes, and its proximity ~ 3.5 Mpc.[] Radio observations at different 
wavelengths |JI| show a structure composed by a compact core, a one- 
sided jet, Double Inner Lobes, a Northern Middle Lobe, and two Giant Outer 
Lobes (See Fig. 10). This morphology, together with the polarization data 
obtained by Junkes et al. |63| and the large-scale radio spectral index dis- 



tribution HI, strongly support the picture of an active radio galaxy with 



a jet forming a relatively small angle (« 30°) with the line of sight. The 
jet would be responsible for the formation of the Northern Inner and Mid- 
dle Lobes when interacting with the interstellar and intergalactic medium, 
respectively. The Northern Middle Lobe can be interpreted as a "working 
surface" [0] at the end of the jet, a place where strong shocks are produced 
by plasma collisions, i.e. it can be considered as the hot spot of a galaxy 
with a peculiar orientation. 



Subluminal velocities of ~ 0.15c have been detected in the jet ||66|| . For 
6 « 30°, we have j3 JET ~ 0.3. One can estimate u ~ 0.4 from the radio 
spectral index of the synchrotron emission in the Northern Middle Lobe and 
the observed degree of linear polarization in the same region. The size of 
the hot spot can be directly measured from the large-scale map obtained by 
Junkes et al. |H| with the assumed distance of 3.5 Mpc, giving as a result 
i? HS ~ 1.75 kpc. A value a ~ 0.01 seems to be reasonable for a source with 



the luminosity of Cen A [58 



The equipartition magnetic field can be obtained for the pre-shock region 
from the detailed radio observations by Burns et al. [62|. The field compo- 
nent parallel to the shock will be amplified in the post-shock region by a com- 
pression factor £. In the case of a strong, nonrelativistic shock front, £ — > 4 
01, and then, if B\\ ~ B ± , we have, B « (£ x / 2 + l) 1 / 2 B ± ^ 5 x 10~ 5 G. 



An enhancement of the Bn component in the Northern Middle Lobe can be 



clearly seen in the polarization maps displayed in Ref. |63 



With the above mentioned values in favor of the input parameters of Eq. 
(30[), the maximum energy of the protons injected in the intergalactic space 

6 Centaurus A (Cen A) was suggested as a possible source of EHECRs by Cavallo ]6l| . 
from quite general energetic arguments. 
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Figure 10: Continuous radio emission from Cen A (1420 MHz). Brightness 
temperature runs from .IK to 50K. The background has been substracted 
(data taken from |63|). 



results £ max ~ 10 21 eV flBBfl. 



We can infer the index in the power-law spectrum from multifrequency 
observations of the synchroton radiation produced by the leptonic component 
of the particles accelerated in the source (see the standard formulae, for 
instance, in the book by Pacholczyk Using the radio spectral index 

obtained by Combi and Romero for the hot spot region ||65| , we get 7 = 2.2. 

Spectral modifications that arise from the interaction of extremely high 
energy protons with the MBR can be computed using Eqs. (|TT| ) and (|I~8|) 
with sufficient accuracy in the case of a nearby source like Cen A. Fig. 11 
shows the modification factor r\. It can be appreciated that the spectrum is 
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Figure 11: Modification factors of southern radio-galaxies. The parameters 
of Pictor A were taken from Ref. 



m 



not dramatically modified, as in the case of most distant sources. 

Finally, it should be mentioned that some observational reports seem to 
back up the above-outlined model. 7-rays in the energy range 300 - 3000 
GeV from Cen A have been detected by Grindlay et al. in 1975 with the 
Narrabi optical intensity interferometer [70|. An excess of cosmic radiation 



from the direction of Cen A in the range 10 
Clay et al. 
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10 17 eV has been reported by 
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3.3 PKS 1333-33 

The radio galaxy PKS 1333-33 is made up of a core, two symmetric jets, 



and two extended radio lobes \72\. The core has been identified with the El 
galaxy IC 4296, which has a redshift z=0. 013 | 73]] . The distance to the source 
is 35.2 h~ l Mpc. The large-scale structure of PKS 1333-33 has been studied 
in detail at 1.3, 2, 6, and 20 cm with an angular resolution of 3.2" W%. The 
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jets are slightly bent, presumably as a consequence of the motion of the core 
with respect to the intergalactic medium. The total flux density of the source 
at 20 cm is ~ 14.5 Jy. The integrated radio luminosity is ~ 2.5 x 10 41 erg 
s -1 assuming a spectral index of a = 0.6 and frequency cut-offs at 10 7 and 
10 11 Hz. Actually, the value a = 0.6 is correct only for the extended radio 
lobes. The spectral index steepens in the jets, reaching a = 1, while the core 
has flatter values: a ~ 0.3. 

An interesting feature of PKS 1333-33 is an intense region of synchrotron 
emission localized at the outer edge of the eastern lobe. Again, this region can 
be considered as a "working surface" formed by the deceleration of the jet. 
This interpretation is supported by VLA polarimetric observations, which 
show a change in the field orientation from parallel to perpendicular to the 



jet axis |72| . This change is probably due to the rearrangement of the field 
lines in the post-shock region. The degree of linear polarization in the eastern 
lobe is in the range 20%-40%. The synchrotron parameters estimated for this 
region are: the minimum energy density (e m i n = 2.1 x 10 -13 erg cm -3 ), the 
minimum magnetic field (-B m m = 1-5 x 10~ 6 G), the minimum total pressure 
(-Pmin = 1-1 x 10~ 13 dyne cm" 2 ), and the degree of linear polarization (m = 
20-40%). 

The leptonic component of the CRs will produce synchrotron emission 
with a spectrum given by S u oc v~ a , where a = (7 — l)/2. Since a = 0.6, we 
get 7 = 2.2. We assume that electrons and protons in the source obey the 
same power law energy spectrum oc E^ 1 . 

The degree of linear polarization expected for the synchrotron radiation 
when the magnetic field is homogeneous is 

m ( 7 ) = f^±| « 71% (32) 
37 + 7 

However, the observed degree of polarization has a mean value of 771(7) ~ 
30%. This fact can be explained by the presence of a turbulent component 
B r in the field, in such a way that 

m W = ™o(7) ( 33 ) 



where B stands for the homogeneous field. From Eq. (R2) and (1331) we get 



B T « 1.2 B (34) 
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and consequently u = B 2 / B 2 otal ^0.6. 

The radius of the acceleration region can be directly measured by means 
of a Gaussian fitting from the detailed VLA maps obtained by Killeen et 



al. |[72|| , resulting R^ s pa 2.5 h kpc. The velocity of the jet is not well 
established. If the source is ~ 10 7 yr old, a velocity /3j et ~ 0.2 can be 
estimated from an analysis of the energy budget |74| . As in the case of Cen 
A, a « 10" 2 . 

Taking the above considerations into account (with a typical value for the 
total magnetic field in the spot of 5_5 ~ 10 f75|), we obtain from Eq. ([30]) 
the maximun injection energy for protons, E max pa 6 x 10 20 eV |[76|| . Thus, 
we have the following proton injection spectrum 

Q(E) oc E- 2 - 2 , E < E < 6 x 10 20 eV. (35) 

Using the formalism presented in Sec. 2.2.1 it is straightforward to com- 
pute the main characteristic of the evolved spectrum. In Fig. 11 we have 
plotted the relation obtained for 77 in the case of PKS 1333-33. The energy 
loss by photomeson production creates the expected cut-off, and the resulting 
ultra high energy CRs (protons and neutrons) pile up just below the thresh- 
old energy of photopion production, forming a bump. However, the energies 
at which the cut-off and the bump appear seem to indicate that CRs above 
100 EeV might be expected from the location in the sky of PKS 1333-33 
(I pa 313.7°, 6^27.7°). 

Whether this scenario is the correct one or not should be answered in a 
few years by the Pierre Auger Southern Observatory |77], [78| (fluorescence 



detector plus ground array) as well as by the future eyes of the OWL j79 
that will deeply watch into the CR-sky. 

4 Hadronic interactions in the final frontier 
of energy 

The energy spectrum beyond 10 15 eV needs to be studied indirectly through 
the extensive air showers (EAS) CRs produce deep in the atmosphere. The 
interpretation of the observed cascades relies strongly on the model of the 
shower development used to simulate the transport of particles through the 
atmosphere. The parts of the shower model related to electromagnetic or 
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weak interactions can be calculated with good accuracy. The hadronic inter- 
action, however, is still subject to large uncertainties. It generally depends 
on Monte Carlo simulations which extrapolate phenomenological models to 
energies well beyond those explored at acellerators. 

There is a couple of quite elaborate models (the dual parton model (DPM) 
and the quark gluon string (QGS) model [|8T|) that provide a complete 



phenomenological description of all facets of soft hadronic processes. These 
models, inspired on 1/N expansion of QCD are also supplemented with gen- 
erally accepted theoretical principles like duality, unitarity, Regge behavior 
and parton structure (for technical details see p2| ). At higher energies, how- 
ever, there is evidence of minijet production |83j and correlation between 
multiplicity per event and transverse momentum per particle [J84| , suggest- 
ing that semihard QCD processes become important in high energy hadronic 
interactions. It is precisely the problem of a proper accounting for semi- 
hard processes the major source of uncertainty of extensive air showers event 
generators. 

Two codes of hadronic interactions with similar underlying physical as- 
sumptions and algorithms tailored for efficient operation to the highest cos- 
mic ray energies are SIBYLL and QGS JET [p9j] . In these codes, the low 
Pt interactions are modeling by the exchange of Pomerons (a hypothetical 
particle with well defined properties, whose precise nature in terms of quarks 
and gluons is not yet completely understood). Regge singularities are used 
to determine the momentum distribution functions of the various sets of con- 
stituents, valence and sea quarks. In the interaction, the hadrons exchange 
very soft gluons, simulated by the production of a single pair QCD strings 
and the subsequent fragmentation into colour neutral hadrons. In QGSJET 
these events also involve exchange of multiple pairs of soft strings. 

As mentioned above, the production of small jets is expected to dominate 
interactions in the cm. energy above y/s ~ 40 TeV. The underlying idea 
behind SIBYLL is that the increase in the cross section is driven by the 
production of minijets [p5 |. The probability distribution for obtaining N 
jet pairs in a collision at energy y/s is computed regarding elastic pp or pp 
scattering as a difractive shadow scattering associated with inelastic processes 
[pq| . The algorithms are tuned to reproduce the central and fragmentation 
regions data up to pp collider energies, and with no further adjustments they 
are extrapolated several orders of magnitude. 

In QGSJET the theory is formulated entirely in terms of Pomeron ex- 
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changes. The basic idea is to replace the soft Pomeron by a so-called "semi- 
hard Pomeron", which is defined to be an ordinary soft Pomeron with the 
middle piece replaced by a QCD parton ladder. Thus, minijets will emerge 
as a part of the "semihard Pomeron" , which is itself the controlling mecha- 
nism for the whole interaction. After performing the energy sharing among 
the soft and semihard Pomerons, and also the sharing among the soft and 
hard pieces of the last one; the number of charged particles in the partonic 
cascade is easily obtained generalizing the method of multiple production of 
hadrons as discussed in the QGS model (soft Pomeron showers) [81] . 



Both, SIBYLL and QGSJET describe particle production in hadron- 
nucleus collisions in a quite similar fashion. The high energy projectile 
undergoes a multiple scattering as formulated in Glauber's approach P7| ], 
particle production comes again after the fragmentation of colorless parton- 
parton chains constructed from the quark content of the interacting hadrons. 
In cases with more than one wounded nucleon in the target, the extra strings 
are connected with sea-quarks in the projectile. This ensures that the inelas- 
ticity in hadron-nucleus collisions is not much larger than that corresponding 
to hadron-hadron collisions. A higher inelastic nuclear stopping power yields 
relatively rapid shower developments which are ruled out by p-nucleus data 



The most direct way to analyze the differences between the models is 
to study the characteristic of the secondaries generated under similar condi- 
tions. For each hadronic code we generate sets of 10 5 collisions in order to 
analyze the secondaries produced by SIBYLL and QGSJET in pp and pA 
(A represents a nucleus target of mass number A = 10) at different projec- 
tile energies. Short lived final state particles are forced to decay according 
with algorithms included in the SIBYLL and QGSJET packages. We have 
recorded the total number of secondaries (baryons, mesons, and gammas), N, 
produced as a result of the interactions. In all the considered cases we found 
that the number of secondaries coming from QGSJET collisions is larger 
than the ones corresponding to the SIBYLL case. For each particle type, 
two-dimension iVx log(-E') distributions were generated. Figure 12 shows 
selected energy distributions for secondaries produced in p-p collisions.f] It is 
easily seen that when the algorithms are extrapolated several orders of mag- 
nitude, the differences between the predicted number of secondaries grow up 



7 For details the reader is referred to |8£ 
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Figure 12: Two-dimension distributions ( Log (E ± ) vs. number of secon- 
daries) obtained from p-p scatterings (incident energies of p 100 TeV - 10 
EeV - 100 EeV downwards). In the left hand side we present the results of 
QGSJET while the right hand side corresponds to the ones of SIBYLL. 
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dramatically. 

Proton induced air showers are generated using AIRES+SIBYLL and 
AIRES+QGSJET.Q Primary energies range from 10 14 eV up to 10 20 5 eV. To 
put into evidence as much as possible the differences between the intrinsic 
mechanism of SIBYLL and QGSJET we have always used the same cross 
sections for hadronic collisions, namely, the AIRES cross section. 

All hadronic collisions with projectile energies below 200 GeV are pro- 



cessed with the Hillas Splitting algorithm ||91|| , and the external collision 
package is invoked for all those collisions with energies above the mentioned 
threshold. It is worthwhile mentioning that for ultra-high energy primaries, 
the low energy collisions represent a little fraction (no more than 10% at 
10 20 ' 5 eV) of the total number of inelastic hadronic processes that take place 
during the shower development. It is also important to stress that the depen- 
dence of the shower observables on the hadronic model is primarily related 
to the first interactions which in all the cases are ultra high energy processes 
involving only the external hadronic models. All shower particles with ener- 
gies above the following thresholds were tracked: 500 keV for gammas, 700 
keV for electrons and positrons, 1 MeV for muons, 1.5 MeV for mesons and 
80 MeV for nucleons and nuclei. The particles were injected at the top of 
the atmosphere (100 km.a.s.l) and the ground level was located at sea level. 

We have analyzed in detail the longitudinal development of the showers. 
The number of different kind of particles have been recorded as a function of 
the vertical depth for a number of different observing levels (more than 100). 

The charged multiplicity, essentially electrons and positrons, is used to 
determine the number of particles and the location of the shower maximum 



by means of four-parameter fits to the Gaisser-Hillas function In Fig. 
13, (A max ) is plotted versus the logarithm of the primary energy for both, 
the SIBYLL and QGSJET cases. It shows up clearly that SIBYLL showers 
present higher values for the depth of the maximum, and that the differences 
between the SIBYLL and QGSJET cases increase with the primary energy. 
This is consistent with the fact that SIBYLL produces less secondaries than 
QGSJET -as discussed before- and as a result, there is a delay in the electro- 
magnetic shower development which is strongly correlated with tt° decays. 
The fluctuations, represented by the error bars, decrease monotonously as 



8 AIRES is a realistic air shower simulation system which includes electromagnetic in- 
teractions algorithms [p0| and links to the mentioned SIBYLL and QGSJET models. 
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Figure 13: Simulation results for the average slant depth of maximum, 
(X max ) , for proton induced showers, plotted versus the logarithm of the pri- 
mary energy. The error bars indicate the standard fluctuations (the RMS 
fluctuations of the means are always smaller than the symbols). The squares 
(circles) correspond to SIBYLL (QGSJET). 



long as the energy increases, passing roughly from 95 g/cm 2 at E — 10 14 eV 
to 70 g/cm 2 at E = 10 20 - 5 eV. 

Using the recorded particle data, we have evaluated lateral distributions 
not only at ground altitude but also at predetermined observing levels. In 
the Figs. 14 and 15 we present the distributions corresponding to a subset 
of all the levels considered, taking into account particles whose distances to 
the shower axis are larger than 50 m. 

The high-altitude lateral distributions show important differences be- 
tween SIBYLL and QGSJET; such differences diminish as long as the shower 
front gets closer to the ground level. The behaviour can be explained taking 
into account the differences between the number of SIBYLL and QGSJET 
secondaries previously reported. Due to the fact that SIBYLL produces less 
number of secondaries, they have -in average- more energy and therefore 
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Figure 14: Comparison between the recorded muon lateral distributions dis- 
played by SIBYLL and QGSJET at different atmospheric altitudes. Incident 
primary energy 10 19 eV. 
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the number of generations of particles undergoing hadronic collisions is in- 
creased with respect to the QGSJET case. As a result, during the shower 
development SIBYLL is called more times than QGSJET, and this generates 
a compensation that tends to reduce the difference in the final number of 
hadronic secondaries produced during the entire shower, and consequently in 
the final decay products, that is, electrons, gammas and muons. 

The lateral distributions of electromagnetic particles are remarkably sim- 
ilar at both (X max ) and ground level.[] However, it comes out from a more 
detailed analysis of the ground distributions that they are not strictly coin- 
cident and that the ratio between SIBYLL and QGSJET predictions does 
depend on r, the distance from the core. In fact, for electrons, this ratio runs 
from 1.25 for small r to 0.73 for r ~ 1000 m, being equal to 1 at r ~ 350 m. 
A similar behaviour is observed for gammas where the lateral distributions 
intersect at r ~ 1000 m. In the case of lateral muon distributions, QGSJET 
predicts a higher density for all distances, but the SIBYLL/QGSJET ratio 
is not constant, ranging from 0.74 near the core to 0.56 at 1000 m. 

Summing up, the most oustanding difference between SIBYLL and QGSJET 
is reflected in the predicted number of secondaries after single p-p and p-nuclei 
collisions. Such a difference increase steeply with rising energy. The different 
number of secondaries predicted remains noticeable during the first stages 
of the shower development, however, the evolution of lateral distributions 
along the longitudinal shower path allows us to clearly observe how the dif- 
ferences in the distributions become monotonously damped, yielding rather 
similar shapes when reaching the ground. Further, we have shown that the 
differences observed at ground level do depend on the distance to the shower 
core. Consequently, we are convinced that it will be possible to obtain rel- 
evant information about the hadronic interactions in air showers from the 
measurement of particle densities at distances far from (as well as close to) 
the shower core. This can be achieved if CR experiments are designed with 
appropiate dynamic ranges. 

Measurement of particle numbers at high atmospheric altitudes (fluorence 
detectors) together with shower maximun and changes rates, should con- 
tribute positively in the understanding of hadronic interactions if the kind 
of cosmic particle that induces the shower could be determined separately. 
Finally, it should be noted that most of the model discrepancies in hadronic 

9 We want to stress that (X max ) is different in each model. 
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Figure 15: Same as Fig. 14 for the case of gammas. 
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interactions will be naturally reduced with the help of data obtained from 
future accelerator experiments like the well known Large Hadron Collider 
(LHC). 

5 Looking Forward 

Cosmic rays of extremely high energy are a weird phenomenon of nature 
for which no truly satisfactory explanation has been found yet. This may 
simply reflect our present ignorance of conditions or processes in some highly 
energetic regions of the Universe - or may imply that exotic mechanisms are 
at play. 

After thirty years of careful work by several groups all over the world, we 
are in possession of a tantalizing body of data, sufficient to raise our curiosity 
and wonder, but not to succeed unravelling in the mistery. The coming 
avalanche of high quality CR-observations promises to make the begining of 
the next millenium an extremely exciting period for CR-physics. 
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